/*!  \file  GeBRD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      09/19/2013 11:18:21 AM
 *   \copyright GNU Public License.
 */
#include <assert.h>
#include <iostream>
#include <cstdlib>
#include "Matrix.h"
#include "UnitTester.h"
#include "ImplementationFactory.h"
#include "BiReduction.h"
#include "GeBRD.h"
#include "setThreads.h"
extern "C" {
void dgebrd_	(	
const int* m,
const int* n,
double 	* a,
const int* lda,
double 	* d,
double 	* e,
double 	* tauq,
double 	* taup,
double 	* work,
const int* lwork,
int 		* info 
);
void dorgbr_ (	
const char   *     ,
const int    * 	m,
const int    * 	n,
const int    * 	k,
      double *    a,
const int    *  lda,
const double *  tau,
double       *  work,
const int    *	lwork,
int          *	 info 
);
}
/*  
 DGEBRD reduces a general real M-by-N matrix A to upper or lower
 bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.

 If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

  The matrices Q and P are represented as products of elementary
  reflectors:

  If m >= n,

     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)

  Each H(i) and G(i) has the form:

     H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T

  where tauq and taup are real scalars, and v and u are real vectors;
  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
  tauq is stored in TAUQ(i) and taup in TAUP(i).

  If m < n,

     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)

  Each H(i) and G(i) has the form:

     H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T

  where tauq and taup are real scalars, and v and u are real vectors;
  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
  tauq is stored in TAUQ(i) and taup in TAUP(i).

  The contents of A on exit are illustrated by the following examples:

  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):

    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
    (  v1  v2  v3  v4  v5 )

  where d and e denote diagonal and off-diagonal elements of B, vi
  denotes an element of the vector defining H(i), and ui an element of
  the vector defining G(i).}
*/
void GeBRD::compute(Matrix & A, Matrix & Q, Matrix & B, Matrix & P ) {
	assert(A.isGeneral()    );
	assert(P.isReflector()  || P.isOrthogonal());
	assert(Q.isReflector()  || Q.isOrthogonal());
	
	assert(B.isBiDiagonal() && B.sameDim(A));
	assert(A.rows() == Q.rows() && A.cols() == P.cols());
	// the thing is that GeBRD actually stored U and V in A 
	//
	if (A.isColMajor()) {
		reset();
		int m   = A.rows();
		int n   = A.cols();
		int min = m < n ? m : n;
		int lda = A.lda();
		double * a = A.getModData();
		double * d = B.getModData();
		double * e = B.getModOffDiagonal(); 
		
		double * taup = P.isReflector() ? P.getModData() : NULL;
		double * tauq = Q.isReflector() ? Q.getModData() : NULL;
	
		if (taup == NULL) allocate(taup, min);
		if (tauq == NULL) allocate(tauq, min);

		double wkopt = 0.0;
		double * work = NULL;
		int lwork = -1, info = 0;
		// query optimal work size
		dgebrd_	(&m, &n, a, &lda, d, e, tauq, taup, &wkopt, &lwork, &info);
		lwork = (int) wkopt;	
		allocate(work, lwork);
		success &= info == 0 && work != NULL;
		if (info == 0 ) {
			if (work) {
				dgebrd_	(&m, &n, a, &lda, d, e, tauq, taup, work, &lwork, &info);
				success &= info == 0;
			}
			
			if ( P.isOrthogonal() ) {
				MISSING;
				char ch = 'P'; 	int reqLWord = -1;
				int m1 = m, n1 = n, k = n;
				P.fill(0);
				cout << "P before " << P <<  endl;
				dorgbr_(&ch, &m1, &n1, &k, a, &lda, taup, &wkopt, &reqLWord, &info);
				cout << "info = " << info << endl;
				reqLWord = (int) wkopt;	
				if (reqLWord > lwork) { // we need more space
					delocate(work);	
					allocate(work, reqLWord);
					lwork = reqLWord;
				}
				dorgbr_(&ch, &m1, &n1, &k, a, &lda, taup, work, &reqLWord, &info);
				cout << "info = " << info << endl;
				delocate(taup);
			} else {
				P.setRefl(A, false);
				P.makeRowMajor();
			}
			if ( Q.isOrthogonal() ) {
				MISSING;
				delocate(tauq);
			} else {
				Q.setRefl(A, true);
				Q.makeRowMajor();
			}
			delocate(work);	
		}
		if (m >= n) B.makeUpper(); else B.makeLower();
	} else {
		A.transpose();
		B.transpose();
		compute(A, P, B, Q);
		B.transpose();
	}
}
GeBRD::GeBRD() : BiReduction("Lapack dgebrd")  {
	setCompileTime(__DATE__, __TIME__);
	canBeParallel = true;
}
bool GeBRD::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
REGISTER(GeBRD);
